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The principle and the efficiency of the Monte Carlo transfer-matrix algorithm are discussed. 
Enhancements of this algorithm are illustrated by applications to several phase transitions in lattice 
spin models. We demonstrate how the statistical noise can be reduced considerably by a similarity 
transformation of the transfer matrix using a variational estimate of its leading eigenvector, in 
analogy with a common practice in various quantum Monte Carlo techniques. Here we take the 
two-dimensional coupled XF-Ising model as an example. Furthermore, we calculate interface free 
energies of finite three-dimensional 0{n) models, for the three cases n — 1, 2 and 3. Application 
of finite-size scaling to the numerical results yields estimates of the critical points of these three 
models. The statistical precision of the estimates is satisfactory for the modest amount of computer 
time spent. 
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, Many important problems in computational physics and chemistry can be reduced to the computation of dominant 
CSj ■ eigenvalues of matrices of high or infinite order. Among the numerous examples of such matrices are quantum 
] mechanical Hamiltonians and transfer matrices. The latter Titere introduced in statistical mechanics by Kramers and 
i: Wannier in 1941 to study the two-dimensional Ising modeU, and ever since, important work on lattice models in 
classical statistical mechanic has been done with transfer matrices, producing both exact and numerical results. 

The analogy of the time-evolution operator in quantum mechanics and the transfer matrix in statistical mechanics 
allows the two fields to share numerous techniques. Specifically, a transfer matrix T of a statistical mechanical 
lattice system in d dimensions often can be interpreted as the evolution operator in discrete, imaginary time t of a 
' quantum mechanical analog, as is well known. That is, T w exp(— tTi), where Ti is the hamiltonian of a system in 
^1 d — 1 dimensions, the quantum mechanical analog of the statistical mechanical system. From this point of view, the 
O ' computation of the partition function and of the ground-state energy are essentially the same problems: finding the 
J-^ ■ largest eigenvalue of T and of exp{—tH), respectively. 

The transfer-matrix Monte Carlo method used in this paper employs an algorithm as simple as the diffusion Monte 
Carlo algorithm, which was developed to compute the dominant eigenvalue of the evolution operator exp(— tTi) . In 
^ ^ contrast to diffusion Monte Carlo, transfer-matrix Monte Carlo provides exact eigenvalues, subject only to statistical 
' noise and as qualified below in Section |l|. More specifically, unlike transfer-matrix Monte Carlo, diffusion Monte Carlo 
suffers from a systematic error, the time-step error, because of the necessity to employ an approximate, short-time 
evolution operator. Sirnilar errors are also found in path-integral Monte Carlo and, in general, in all approaches based 
on the Trotter formulacl. An alternative, related approach, viz. Green function Monte Carlo, used to compute the 
dominant eigenvalue of {Ti — E)~-^, where E is close to the ground state energy, does not suffer from a time-step error, 
and, from that point of view. Green function Monte Carlo is more elegant than diffusion Monte Carlo. However, the 
Green function Monte Carlo algorithm is considerably more complicated, and enhancement of that algorithm by the 
variance reduction techniques discussed below, has its limitations. 

From an orthodox complexity theory point of view, exact numerical transfer-matrix computations for lattices in 
more than one dimension are intractable, since the order of transfer matrices grows exponentially with the number of 
lattice sites in a transfer slice. Standard Monte Carlo methods in statistical mechanics, on the other hand, statistically 
sample the Boltzmann distribution, typically employing some variant of the Metropolis algorithm. One can argue 
that Monte Carlo methods are of polynomial complexity in the system size, at least for certain important physical 
observables. This raises the question of the ultimate utility of the transfer matrix for computational purposes. 

In many cases, one is interested in the behavior of systems in the thermodynamic limit. For critical systems 
in particular, one has to rely on finite-size scaling and extrapolation methods to extract the relevant information 
from the computations. The transfer-matrix method has advantages in both respects. Firstly, one can compute the 
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spectrum of the transfer-matrix method virtually to machine precision, which permits extrapolation without serious 
loss of numerical accuracy. Secondly, a large body of numerical evidence suggests that the transfer-matrix spectrum 
has weaker corrections-to-scaling than quantities commonly computed by standard Monte Carlo. Clearly, also the 
transfer-matrix Monte Carlo method takes advantage of the weakness of the corrections-to-scaling. Unfortunately, 
statistical noise is introduced, but this can be substantially reduced by the use of optimized trial eigenvectors, by 
virtue of which the Monte Carlo process is in effect only used to compute corrections to an already sophisticated 
approximation. 

If one could neglect the correlations introduced by the re-weighting step of the transfer-matrix Monte Carlo algo- 
rithm [see the split /join steps ( pa| ) and (2b) in the algorithm given in Section ||] and if one could ignore the resulting 
loss of efficiency of the transfer-matrix Monte Carlo algorithm, this method would be a solution to the exponential 
growth problem mentioned aboveu. In addition, transfer-matrix Monte Carlo would be completely free of critical 
slowing down, since the correlation time of the algorithm is equal to the correlation length of the slices used in the 
definition of the transfer matrix. Again, the use of optimized trial eigenvectors can serve to reduce the detrimental 
effect of the multiplicative re-weighting. 

Another feature of the Monte Carlo transfer matrix, which can contribute to a reduction of the correlation time of 
the stochastic process, is that moves are effectively made at surface sites. This makes it much easier to overcome the 
barriers some systems present to standard Monte Carlo algorithms. An example of such a system is the XK-Ising 
model discussed in Ref. ^ 

The layout of this paper is as follows. In Section || we review the basic Monte Carlo algorithm to determine 
transfer-matrix eigenvalues by means of a statistical implementation of the power method. Apart fr om relatively 
minor details, the algorithm given in Section || is the same as the one discussed in Refs. ||, ^ 0- Section III describes 
the similarity transformation of the transfer matrix, which leads to a pronounced decrease of the statistical errors of 
the Monte Carlo process. Section III in particular describes in detail the construction of a variational approximation of 
the eigenstate associated with the largest eigenvalue. This approximate eigenstate yields the similarity transformation 
used to reduce the statistical errors of the algorithm. Details of the speed-up of the algorithm are presented at the 
end of Section III a coupled Xy-Ising model in two dimensions. Finally, Section |^ contains applications of the 
transfer-matrix Monte Carlo method to three-dimensional 0{n) models ipjrfi = 1,2 and 3. Preliminary discussions 
of the the work discussed in Sections HI and IV were published elsewhereQ'cl. 



II. MONTE CARLO IMPLEMENTATION OF THE POWER METHOD 



Consider an operator T of which we want to compute the dominant eigenvalue. Let T be represented by matrix 
elements (i?|T|S') — Trs, where \R) and \S) are basis states of the physical system under consideration. These states 
will be treated here as discrete. For Monte Carlo calculations, the distinction between continuous and discrete states 
is a minor technicality; in the discussion below, generalization to the continuous case follows immediately by replacing 
the appropriate sums by integrals and replacing Kronecker by Dirac 5- functions. 

Perhaps the simplest way to calculate the dominant eigenvalue of a matrix or integral kernel is the power method. 
That is, choose an arbitrary initial state and compute iteratively: 

tu(*+i)} = — T|u«), (1) 

Ct+l 

where Ct+i is a constant chosen so that |7i(*+^-') is normalized or in some other convenient standard form. For t oo, 
the constants Ct approximate the dominant eigenvalue Aq of T and the vectors ju^*)) converge to the corresponding 
eigenvector. 

To implement Eq. (|l|) by Monte Carlo, ju^*)) is represented by a sequence of Nt walkers. Each of these walkers is a 
pair (Ra, Wa), a = 1, . . . , Nt- The variable Ra of a walker represents a possible configuration of the system described 
by T, and Wa represents its statistical weight. The latter quantity is subject to the condition wi < Wa < w^, where 
wi and Wu are bounds introduced so as to keep all weights Wa of the same order of magnitude, which improves the 
efhciency of the algorithm. This sequence of walkers represents a (sparse) vector with components 

Nt 

Ur ^^WaSR,B^„, (2) 
a = l 

where S is the usual Kronecker (5-function. The underscore is used to indicate that the j/^ represent a stochastic 
vector lu*-*-*). A stochastic process will be defined presently with transition probabilities such that Q+i has a 
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conditional expectation value equal to T|m*^*)) for any given sequence of walkers representing !«*•*•'). In practice, one 
has to average over the stationary state of a stochastic process in which the constants Ct are determined on the fly, so 
that Ct+i and are correlated. As a consequence, there is no guarantee that the stationary state expectation 

value of is precisely an eigenstate of T, at least not for finite Nt- The same mathematical pttpblem occurs if 

one takes the time-average of Eq. (0) in the presejnce of noise correlated to the Ct ■ The resulting biasBB has also been 
discussed in the context of diffusion Monte CarloQ. 

To define the stochastic process, Eq. dll) is rewritten as 



(t+i) 

^ Ct+l 



= —Y.PrsDsu^'^S, (3) 



where 



Ds^Y. '^RS and Prs = Trs/Ds- (4) 

R 

Eq. (|^) describes a process represented by a Monte Carlo run which, in addition to a few initial equilibration sweeps, 
consists of a time series of a little over Mq sweeps over all walkers at times labeled hj t — . . . , 0, 1, . . . , Mq. The sweep 
at time t consists of tjwo steps designed to perform stochastically the matrix multiplications in Eq. (|^). Following 
Nightingale and Blotctl, the process is defined by the following steps, which transform the generation of walkers at 
time t into the the generation at time i + 1. Variables pertaining to times t and t + 1 will be denoted respectively by 
unprimed and primed symbols. 

1. Update the old walker {Sa,Wa) to yield a temporary walker {S'^,w'^) according to the transition probability 
Ps'^Sai where w'^ — Ds^Wa/c' , for a — I, ... , Nt- The next step can change the number of walkers. To maintain 
their number close to a target number, say A'o, choose c' — Xo{Nt/NQy^'^ , where Ao is a running estimate of the 
eigenvalue Aq to be calculated, where s > 1 (see below). 

2. From the temporary walkers construct the new generation of walkers as follows: 

(a) Split each walker {S'jw') for which w' > into two walkers (S", ^w'). The choice &u = 2 is a reasonable 
one. 

(b) Join pairs and {R'p,w'p) with w'^ < bi and < bi to produce a single walker {R'^,w'^ + w'^), 
where Ri.^ — R'^ or R'^ — Sp with relative probabilities w'^ and w'^. We chose b\ = 1/2. 



(c) Any temporary walker left single in step (2b), or for which bi < w'^ < b^, becomes a permanent member of 
the new generation of walkers. 

The algorithm described above was constructed so that for any given realization of | u^*) ) , the expectation value of 
Ct+i\u^*^^^), in accordance with Eq. ([|), satisfies 

e(q+i|^.(*+i))) =T|jiW), (5) 

where E(-) denotes the conditional average over the transitions defined by the above stochastic process. More generally 
by p-fold iteration one findsEl: 



E 



p 

.6=1 



u 



(*+p)) =TP|uW}, (6) 



The stationary state average of is close to the dominant eigenvector of T, but, as mentioned above, it has a 
systematic bias when the number Nt of walkers is finite. For increasing p, components of non-dominant eigenvectors 
can be projected out and thus the bias is reduced, in principle. Unfortunately, the variance of the corresponding 
estimators increases as their bias decreases. The reader is referred to Refs. ||, |[, ^ for a more detailed discussions 
of this problem. Suffice it to mention here, firstly, that s is the expected number of time steps it takes to restore the 
number of walkers to its target value A^orapd, secondly, that strong population control (s = 1) tends to introduce a 
stronger bias than weaker control (s > 

With Eq. (||) one constructs an estimatorO of the dominant eigenvector of the matrix T: 



Mo /p-l \ 
^ t=l \6=0 / 



)■ (7) 
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More practically, suppose that {tpT] is an approximate leading eigenbra (■0t| of T, and that O is an arbitrary 
operator. The mixed expectation value of O can be approximated as 

An important special case is obtained by choosing in this expression O = T and {iPt\R) — 1 for all R. The latter 
corresponds to the infinite-temperature approximation for the trial state and in that case, Eq. (|^) reduces to an 
estimator for the dominant eigenvalue of T: 



(nCo ^ 



where 



Wt = (V'tIm^*^) = (10) 

For the above special choice of the trial bra {ipT \ , Eq. (||) becomes the expression for the surface expectation value 
of O in the geometry shown on the right in Fig. ^ Although we have used the transfer-matrix algorithm only for the 
computation of the dominant eigenvalue of the transfer matrix for the applications discussed in this paper, it should 
be mentioned for completeness that one can also compute bulk expectation values, at least asymptotically, as follows. 

One can represent the Kramers- Wannier transfer matrix by the graph shown in Fig. ^.a. This matrix transfers from 
an old slice to a new one, with slices represented respectively by small full and large open circles. The process adds 
only one new site: the open circle labeled 1. One site, the small closed circle labeled L, is about to disappear into the 
bulk. Coincidences of both types of circles represent Kronecker-(5 functions in the transfer matrix [see Eq. j]!^)]. The 
solid lines stand for interactions added in one transfer operation. One can define a transfer matrix with extended slices 
consisting of m of the original, minimal slices. The dominant eigenvector of this extended transfer matrix is simply 
the original eigenvector multiplied by the Boltzmann weight associated with the portion of the lattice containing 
variables that have not yet been summed over. Eq. (^, used with any operator in which occur only variables of slice 
TO, becomes a bulk expectation value .for to — > oo. The implementation of this concept is called forward walking in the 
context of quantum Monte CarloMtJ, and this only requires extending the walkers so that their states correspond to 
the extended slices introduced above. This increases the memory requirements and the cost of splitting a walker, but 
otherwise the efficiency of the algorithm is not affected. 



III. VARIANCE REDUCTION (IMPORTANCE SAMPLING) AND TRIAL VECTORS 

In principle, if {iPt\ equals an exact eigenbra of the operator O in equation Eq. (H), the right-hand side of the 
expression is a zero- variance estimator. In general, no exact eigenvectors are known, but even an approximation may 
yield a substantial reduction of statistical noise. A more efficient well-knownEj way to exploit an approximate left 
eigenbra {tpT \ to reduce variance works by application of the method described above to a similarity transform of the 
original operator T. This transformation is defined by: 

T = ITI"\ (11) 

where I is diagonal in the configuration presentation, and is defined as 

I = J2\R){i:T\R){R\. (12) 

R 

Ideally, (V'tI would equal the exact dominant eigenbra of T. In that case, the stochastic process defined as above, 
but with T replaced by T, would become optimally efficient and in fact would lack critical slowing down. For such 
an ideal process D, defined as in Eq. (jj) as a function of T, would be a constant times the unit matrix. The walker 
weights would no longer fluctuate so that birth and death process would no longer occur. The walkers would evolve 
into a statistically independent ensemble. The estimator given in Eq. (^), appropriately transformed, would have 
zero variance. The transformed bra {i/jt\ = (V't|I~^ would have all elements equal to unity in the configuration 
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representation. In other words, T would be represented by a stochastic matrix, which would eliminate re- weighting 
of walkers and the concomitant split /join step in the algorithm. 

In the absence of exact eigenbras, approximations may be obtained by variational methods. The variational ex- 
pression for the leading eigenbra {tpT \ can conveniently be cast in the form of an effective surface Hamiltonian with 
pair interactions between nearest neighbors, next-nearest neighbors, and so on. These-Lnteractions are treated as 
variational parameters and can be determined from an analysis of the walker populationliJ. 

Since generalization to higher dimensions and models with different microscopic variables is straightforward, it 
will suffice to consider the Kramers- Wannier transfer matrix for the two-dimensional Ising model to explain the 
construction of trial vectors used in the applications discussed in Section 

For a simple quadratic lattice of M sites, wrapped on a cylinder with a circumference of L spins and helical boundary 
conditions, the transfer matrix for the Ising model is 

T5.K = c^(^i'^i+^^'^-)[]<5,.,,^+,, (13) 

with S = (si, S2, . . . , sl) and R — (ri, r2, . . . , tl), where the Si = ±1 and ri ± 1. The conditional partition function 
of the lattice of M sites, subject to the restriction that the spins on the left-hand edge are in state R, as illustrated 
in Figure |l|, is denoted Zm (R) ■ One has 

ZM+iiS) = J2 Ts,rZm{R). (14) 

R 

Obviously, for M — > oo the restricted sums Zm{R) are proportional to the components u^'' of the dominant right 
eigenvector of the transfer matrix. The eigenvector is represented by the graph on the right in Figure l]. Full circles 
indicate spins that have been summed over, while the fixed surface spins are represented by the open circles; each bond 
represents a factor e:icp{KsiSj). The left eigenvector, which is the one that has to be approximated by an optimized 
trial vector, is represented by the graph on the left. In passing, we mention the following relation between left and 
right eigenvectors, which follows by inspection of the graphs: 



L-l 



\S) = W e^"'*'+i(C/(5)|M(°°)), (15) 



where U is the reflection operator: U{S) = {sl, sl^i, ■ ■ ■ , si). 

A similarity transformation of the transfer matrix T can be introduced by dividing up the interaction energies 
between the columns differently. That is, h is introduced by writing 

T5,fl^e''(^'«). (16) 

A transformation h ^ h is defined by 

h{S,R) ^ g{S) + h{S,R) ~ g{R), (17) 

fs,R^i^T{S)Ts,R/MR), (18) 
^T(5)=es(^). (19) 

For purposes of variance reduction, versatile trial vectors that capture some of the essential physics without seriously 
slowing down computations, can be chosen of the following form 

^^(S) = e^l^^""'"\ (20) 

a form reminiscent of the Jastrow functions used for quantum many-body systems. The asterisk in the sum over pairs 
indicates that the Kij are truncated for distances greater than a couple of lattice spacings. 

The couplings Kij in Eq. (|2^) are variational parameters. They can be determined efficiently with the Monte 
Carlo scheme introduced by Umrigar, Wilson and Wilkinso, i.e., by minimization of the variance of D{S), where the 
variance is approximated by a weighted sum over the states of the walkers of one generation, during the initial stage 
of the Monte Carlo run. This procedure is efficient and stable as long as the Kij are truncated with care, in which 
case it is perfectly feasible to use as many as 50 to 100 different parameters. 



5 



o 



<> 



(a) 



2 
1 



(b) 



FIG. 1. Illustration of left and right eigenvectors of the transfer matrix. 
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FIG. 2. Illustration of the calculation of correlation functions involving spins in the bulk below the surface layer. Site labels before the 
addition of the new spin (open circle) appear to the right, and the new labels to the left of a lattice point. 
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The magnitude of the Kij is expected to increase with the strength of the correlations between surface spins. Since 
all correlations between surface spins for the left eigenvector have to be propagated through the lattice on the left, as 
illustrated in Figure |l|, one expects that for high temperatures, i.e., small K, 

K,,ozK''^^, (21) 

where dij is the length of the shortest path along edges connected by bonds between sites i and j. By inspection 
of the graph in Figure |^, we therefore expect the following partial ordering in decreasing strength of interaction and 
increasing dij 



diL = 2, 



di2 


= C?23 = • 


.^dL 


-1,L 


= dix 


-1 


= 3, 


di3 


= C?24 = ■ 


■ = dL 




= dix 


-2 


= 4 


dii 


= C?25 = ■ 


■ = dL 


-3,L 


= dix 


-3 


= d2,L = 5 



(22) 



It is important to note that if Kij = A'^+ij+i the corresponding factors cancel in the transformed transfer matrix 
T for 2 < i < L — 2, since Si = i^+i for non-vanishing transfer-matrix elements. For reasons of efficiency it is therefore 
advantageous to have this equality satisfied as often as possible. Unfortunately, helical boundary conditions introduce 
a step which destroys translation symmetry on the surface and renders the partial ordering in Eq. ( p^ ) insufficient. 
For example, sites 1 and 2 are more strongly correlated than 2 and 3, and correlations keep decreasing through pair 
(i - 1, L). Consequently, K12 > K23 > . .. > Kl^i^l- 

In practice, the differences between the Kij with dij = 3 are frequently greater than the higher-order Kij. Then, it 
is necessary to treat K12 and K23 as different parameters of the trial vector. An efficient compromise is to treat Kij 
in which site 1 or L participate as different. The same applies to all Kij for which the shortest path between i and 
j straddles the step on the surface. To summarize, we distinguish different types of pairs of sites both on the 
basis of the distance dij and to some extent on the location of the pair, enforcing as much translation invariance as 
possible. 

Clearly, none of the above depends only on lattice geometry or the Ising nature of the variables. In general, only a 
method is required to generate lists of lattice sites separated by various distances dij . These can be constructed with 
simple graph theoretic tools such as incidence matrices, which makes it possible to deal with different dimensions and 
lattice types in an identical fashion once the pertinent incidence matrix has been defined. 

To illustrate the efficiency and flexibility of this technique for constructing trial vectors, we use the XK-Ising model. 
It consists of coupled Ising and planar rotator degrees of freedom on a simple quadratic lattice. On each lattice site 
there are two variables: Si ± 1 and rii, a two-component unit vector. The reduced Hamiltonian — divided by —k^T — 
is given by 

H — {A Mi ■ rij + B rii ■ rijSiSj + CsiSj) . (23) 

We consider the special case A — B and only from the point of view of the performance of transfer-matrix Monte 
Carlo algorithm. For a discussion of the physics of this model the reader is referred to Ref. |^. The trial vectors 
discussed above for the Ising model have an immediate generalization: 

ipT = exp {Aij Ui ■ rij -I- B^j • rijSiS-j -t- Ci^SiSj) | . (24) 

The truncation scheme introduced above for the Ising model is purely geometrical, and therefore carries over without 
changes to the XY-lsing model. It should, however, be noted that there are models and choices of transfer matrices 
to which the above scheme is not applicable. Ref. ^ contains a discussion and an example of such a case. 

Table | shows the estimates of the dominant eigenvalue of XY-ls'mg model for trial vector truncated at different 
values of dij. As can be seen by comparing the first and last lines of the table, the variance in the estimate of the 
eigenvalue is reduced by a factor 300 for a fixed number of Monte Carlo steps. Taking into account that the computer 
time per step doubles, this constitutes a speed-up by a factor of 150. 
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TABLE I. Estimated eigenvalue and standard deviations for the JfY-Ising model. These data apply to the point 
(A = 1.005, C = -0.2285) [c/. Eq. @] on the line where Ising and XY transitions coincide. Results are shown for various values 
of dm, the path length of the cutoff in Eq. (p^). The results are for a strip of width L = 20 and were obtained with a target number of 
walkers Nq = 10, 000 and Mq = 1, 250Z; generations of which an initial 10% were discarded. The last column shows the computer time in 
arbitrary units needed per time step of one walker. 
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IV. APPLICATIONS 

As an illustration of the transfer-matrix technique we apply the method to three-dimensional 0{n) models for 
n — 1, 2 and 3, i.e. the Ising, planar and Heisenberg model. In particular the significance of the results for of the 
planar and Heisenberg models goes beyond mere illustrations. These results are sufficiently accurate to be of some 
relevance for the location of the critical points. 

The 0(n) spins are located on the simple cubic lattice. The transfer matrix for an L x L x cx) system, with helical 
boundary conditions and layers of N = sites each, is a straightforward generalization of Eq. ([l3| ) and reads: 

N-l 

Ts,R = n '^^-"■,+1 exp [Ksi ■ (n + tl + rjv)] , (25) 

1=1 

where the and are n-component unit vectors, S, — (si, S2, . . . , s^r) and R, = (ri, r2, . . . , r^r). 

As discussed above, transfer-matrix Monte Carlo is designed to compute the dominant eigenvalue Aq of the transfer 
matrix. The reduced free energy per site is / = — In Aq . From the free energy one can calculate the surface tension as 
the difference in free energy of two systems: one with ferromagnetic and the other with antiferromagnetic interactions, 
if the dimensions are chosen so as to force an interface in the antiferromagnetic system. For L x L x oo systems with 
helical boundary conditions, to which the present calculations are restricted, this means that L has to be even. 

Renormalization group theory predicts that the values of A, the reduced interface free energy per lattice site, as a 
function of coupling K and system sizes L collapse onto a single curve, at least close to the critical point and for 
sufficiently big systems. In terms of the non-linear thermal scaling field 

u{K) = K -K^ + a{K - K^f + ..., (26) 

this curve S(x) is determined by 

l\{u,L) = L^-'^Y.{Ly^u), (27) 
for a d-dimensional system with a thermal scaling exponent j/t- The function E can be expanded in a series: 

oo 

S(a;) = ^E/a;', (28) 
;=o 

and for 0{n) models behaves for large x as: 

E(a;) = (29) 

where p(_l) = ^'^'^ P(2) = p(3) — 1. 

Eqs. ( p?! ) to (^) are useful for the interpretation of the 0(n) transfer- matrix Monte Carlo results for the interface 
free energy. These results were obtained using finite sizes up to L = 12, and populations typically consisting of 2500 
or 5000 walkers. Typical run lengths are 5000 steps, where each step means the addition of a surface layer of L x L 
spins. Variance-reducing trial vectors [see Eq. (|2^)] were constructed for path lengths up to 5. As before, the variance 
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FIG. 3. Finite-size scaling plot for the interface free energy of the three dimensional Ising model. 

of the Monte Carlo process was observed to decrease considerably with increasing path length. For each system size, 
interface free energies were obtained for approximately 10 different couplings in a range of about 10% around the 
critical points of the Ising and planar models, and about 1% for the case of the Heisenberg model. 

On the basis of these results for the Ising (n = 1) case, the function E is shown in Figure ^ This data collapse is 
achieved by means of a least-S6||Uares fit with parameters Kc, ut, a, and 13 Taylor coefficients S/, a generalization 
of a technique used in the pasiu. 

To check if the system sizes were in the asymptotic finite-size scaling regime, fits were done both with and without 
the 6 X 6 X cx) data. The results of these fits are displayed in Table |ll| in Appendix A. To summarize, the results 
are: = 0.22162 ± 0.00002 and yr = 1.584 ± 0.004 using data with i = 6 through 12; and = 0.22167 ± 0.00004 
and t/T = 1.584 ± 0.014 if the L — 6 are omitted. These results agree well with accurate determinations using other 
methods (see e.g. Ref. p^and references therein) which appear to cluster about Kc = 0.221655 (with a margin 

of about 10~^) and yr — 1.586 (with a precision of a few times 10~^). 

It is remarkable that the corrections to scaling appear to be very small, as appears from the data shown in Figure 
In standard Monte Carlo analysest^l of L x L x L systems these corrections are quite prominent, and form an obstacle 
to the accurate determination of critical parameters. 

The scaling plot shown in Figure ^, can be used to determine the amplitude As graphically: on a double logarithmic 
plot the asymptotic slope of the curve follows from the known value of the thermal exponent j/t, cf. Eq. (p9|). The 
problem of calculating this-, amplitude has attracted considerable attention lately and the reader is referred to a 
paper by Shaw and FishcrEj for details and further references to the literature. For the largest values of the scaled 
temperature variable x, we fiad A'^^ ~ A-^Kc^^'^ — 1.8, while the trend with x is an increasing one. This value is 
somewhat larger than Mon'sEll estimate A-^ — 1.58 ± 0.05, but still in the range 1.4 < A-^ < 2.0 obtained by Shaw 
and Fisher. As a final comment we note that Mon's method requires systems of linear dimensions in excess of 48 to 
reach the asymptotic infinite-size regime, with an increasing trend of the estimates of As with increasing x = L^^u. 
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FIG. 4. Finite-size scaling plot for the interface free energy of the three dimensional planar model. 



A similar analysis was performed for the planar model (n = 2). In comparison with the Ising case, the scaling 
function E behaves more smoothly as a function of x, so that a satisfactory fit could be obtained with fewer Taylor 
coefficients. The fitted parameters, which are K^, j/t, a, and 8 Taylor coefficients I]„, are shown in Table III 
of Appendix A. Our results for the critical point are Kc — 0.45410 ± 0.00003 for system sizeSpLj-^ 6 to 12, and 
Kc = 0.45413 ± 0.00005 for i = 8 to 12,-|These values are close to results from series expansion£3'El A'c = 0.45386 
and standard Monte Carlo calculationaHJ Kc = 0.4531 (no errors quoted). Also our results for the temperature 
exponent, namely j/t = 1.491 ± 0.003 for L > 6 and j/t = 1-4^ ± 0.006 for L > 8 are in a good agreement with 
existing results; we quote the coupling-constant-expansion valueO j/t = 1.495 ± 0.005. 

Fitted with these parameters the data collapse very well onto the function S, as shown in Figure ^. Again, this 
scaling plot can be used to determine the amplitude As graphically: in this case the asymptotic power-law exponent 
is l/j/T- A fit of the data at the highest available values of x = L^^u leads to Aj^ — 5.9, while the trend is still 
increasing with x. 

The calculations for the Heisenberg case n = 3 were clustered in a narrow interval around the critical temperature, 
and were not aimed at an accurate determination y^. Thus, the transfer-matrix Monte Carlo data could be analyzed 
by means of a least-square fit with less parameters: K^, j/t and 3 Taylor coefficients E„. The fit is shown in Table IV in 
Appoiidix A. The result for j/t is well within the statistical accuracy, equal to the known coupling-constant-expansion 
valueEa j/T = 1.418. Including the latter value as a known variable in the fits leaves our results for the critical point 
practically unchanged. These are: = 0.69291 ± 0.00004 for system sizes L-f^ 6 to 12, and = 0.69294 ± 0.00008 
for L = 8 to 12. These values are close to results from series expansion^J: Kc — 0.6916, and more recentlynJ: 
Kc = 0.69294; and from Monte Carlo calculationsEZi: Kc = 0.693035 ± 0.000037. The difference with our result with 
the L — 6 data included could be interpreted as an indication of a small finite-size effect. 

The data collapse for the n = 3 case onto the function E as determined by the least-squares fit is shown in Figure |^. 

Finally we remark that, although in each of the cases n = 1, 2 and 3 the finite-size effect appears to be small for 
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FIG. 5. Finite-size scaling plot for the interface free energy of the three dimensional Heisenberg model. 
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TABLE II. Parameters, as defined in tlie text, and tiieir standard errors for the scaling function of the interfacial free energy of the 
three-dimensional Ising model. 





n > 6 


n > 8 




0.22162 ± 0.00002 


0.22165 ± 0.00003 


VT 


1.583 ± 0.004 


1.594 ± 0.009 


So 


0.6171 ± 0.0007 


0.6194 ± 0.0025 


El 


2.6111 ± 0.0176 


2.5650 ± 0.0505 


Sa 


6.0475 ± 0.1001 


5.8047 ± 0.2565 


S3 


9.2362 ± 0.3052 


8.4073 ± 0.6724 


S4 


6.0087 ± 0.6350 


4.8601 ± 1.0249 


Ss 


-13.6165 ± 0.9830 


-10.8331 ± 1.7414 


Se 


-33.2578 ± 4.0574 


-24.7108 ± 5.6481 


S7 


4.4790 ± 3.3849 


1.4272 ± 4.2080 


Sg 


70.6918 ± 11.4018 


46.6641 ± 14.0107 


S9 


28.0314 ± 7.9088 


21.4636 ± 8.5327 


Sio 


-69.1548 ± 14.3387 


-40.6066 ± 15.7549 


Sii 


-43.4754 ± 10.0932 


-27.7038 ± 10.3867 


S12 


25.2137 ± 6.4840 


13.2245 ± 6.4356 


Sl3 


18.8658 ± 4.8588 


10.5010 ± 4.6564 


a 


-2.65 ± 0.16 


-2.65 ± 0.37 



TABLE III. Parameters, as defined in the text, and their standard errors for the scaling function of the interfacial free energy of the 
three-dimensional planar model. 





n > 6 


n > 8 




0.45410 ± 0.00003 


0.45413 ± 0.00005 


yr 


1.491 ± 0.003 


1.487 ± 0.006 


So 


1.2448 ± 0.0010 


1.2469 ± 0.0033 


Si 


2.5592 ± 0.0144 


2.5929 ± 0.0345 


S2 


2.4285 ± 0.0439 


2.4738 ± 0.0796 


S3 


0.9881 ± 0.0623 


0.9544 ± 0.1031 


S4 


-0.5096 ± 0.0664 


-0.4292 ± 0.1050 


S5 


-0.7770 ± 0.1579 


-0.5171 ± 0.2741 


Se 


0.1737 ± 0.0754 


0.0793 ± 0.1204 


S7 


0.4346 ± 0.1503 


0.1847 ± 0.2567 


a 


-0.7805 ± 0.1159 


-0.9245 ± 0.2259 



L > 6, it is large for i = 4. For this reason the L = 4 data were not included in the fits. 
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APPENDIX A: SCALING PLOT PARAMETER ESTIMATES 



Tables || through IV contain estimates of the parameters used in the finite-size scaling plots for the interface free 
energy of 0(n) models, as discussed in Section IV. 
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TABLE IV. Parameters, as defined in tiie text, and their standard errors for tiie scaling function of the interfacial free energy of the 
three-dimensional Heisenberg model. The Monte Carlo data were taken relatively close to K^, so that the temperature exponent yx is not 
accurately determined. The accuracy of Kc is unaffected. 





n > 6 


n > 8 




0.69291 ± 0.00004 


0.69294 ± 0.00008 


y-T 


1.44 ± 0.07 


1.55 ± 0.18 


So 


1.8919 ± 0.0015 


1.8933 ± 0.0043 


Si 


2.4563 ± 0.3123 


1.9036 ± 0.7665 


E2 


0.7991 ± 0.3005 


0.5097 ± 0.4651 
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